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The properties of a balanced two-component Fermi gas in a one-dimensional harmonic trap are 
studied by means of the coupled cluster method. For few fermions we recover the results of exact 
diagonalization, yet with this method we are able to study much larger systems. We compute 
the energy, the chemical potential, the pairing gap, and the density profile of the trapped clouds, 
smoothly mapping the crossover between the few-body and many-body limits. The energy is found 
to converge surprisingly rapidly to the many-body result for every value of the interaction strength. 
Many more particles are instead needed to give rise to the non-analytic behavior of the pairing gap, 
and to smoothen the pronounced even-odd oscillations of the chemical potential induced by the shell 
structure of the trap. 


Ultracold gases are ideal systems for engineering highly 
non-trivial states of matter. They allow one to prepare, 
manipulate, and measure with great accuracy strongly 
correlated quantum systems [1-3]. Thus, they pro¬ 
vide a perfect playground for classical simulations, and 
can accurately serve as quantum simulators [4]. One¬ 
dimensional (ID) systems are particularly fascinating be¬ 
cause of the important role played by quantum fluctua¬ 
tions [5]. Experimental studies of strongly correlated ul¬ 
tracold atomic gases started with the seminal works on 
the Tonks-Girardeau gas [6-8]. More recently the su¬ 
per Tonks-Girardeau regime has been achieved [9], and 
the hrst experiments with fermions with a tunable spin 
have been launched [10]. For this Rapid Communication, 
particularly relevant are experiments on two-component 
Fermi gas performed in ID harmonic traps containing 
very few atoms, where the number of spin-up and spin- 
down atoms may be fully, and separately, controlled [11- 
15]. The theoretical studies of ID Fermi gases in the 
many-body regime have a long tradition [16-27]. More 
recently, several papers have addressed in detail the few- 
body case, and its evolution towards the many-body 
regime [28-37]. 

In order to address increasingly more elaborate ex¬ 
perimental findings, the efficient numerical treatment of 
many-body quantum systems stands as one of the great 
challenges of modern physics. Despite enormous progress 
and the development of many powerful approaches (cf., 
e.g., density functional theory [38], exact diagonalization 
[39, 40], quantum Monte Carlo (QMC) [26, 40, 41], den¬ 
sity matrix renormalization group (DMRG) and tensor 
network states [42, 43]), new numerical approaches that 
are able to investigate hitherto unexplored phenomena 
are always more than welcome. 

In this Rapid Communication, we study a balanced 
two-component Fermi gas in a one-dimensional harmonic 


trap by means of a quantum chemistry approach—the 
coupled cluster (CC) method [44-51]. In condensed mat¬ 
ter, this method has up to now been successfully ap¬ 
plied to spin-I/2 lattice models in ID and 2D (see Refs. 
[52, 53] and references therein), and to trapped ultra¬ 
cold bosons [54, 55]. However, the CC method is ideally 
suited to study fermionic systems, where the number of 
occupied orbitals grows at least as fast as the number of 
particles in the system, even in the absence of interac¬ 
tions. CC recovers results known from exact diagonal¬ 
ization [35] and a path integral approach [56] , but it also 
allows one to study much larger systems (up to ~ 80 par¬ 
ticles). QMC methods permit one to look at even bigger 
clouds [26]; ground state properties may be studied by 
means of the very accurate diffusion QMC, while finite 
temperatures may be considered via path-integral MC. 
Finally, CC compares well with the state-of-art DMRG 
calculations [27], with the advantage however of being 
a method explicitly working in continuous space rather 
than in a lattice. 

Here we compute with high precision the energy, the 
chemical potential, the pairing gap and the density pro¬ 
file of the fermionic gas as a function of the number 
of particles and of the interaction strength. While the 
ground state energy of the system converges astonish¬ 
ingly rapidly to the many-body result for any interaction 
strength, the other quantities depend more sensitively 
on system size and interaction strength. In particular, 
strong even/odd oscillations in the chemical potential of 
the trapped gas persist up to a very large number of par¬ 
ticles, and the non-analytic behavior of the BCS pairing 
gap emerges only relatively slowly with the system size. 

Model. We consider a two-component Fermi gas con¬ 
taining N = -I- Ni atoms in a balanced configuration, 

i.e., with = N^. All atoms have the same mass m, 
and they are bound to move along one dimension due to 
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Figure 1. Rescaled energies 8m = Em/Em^ for various num¬ 
ber of fermions, as a function of the interaction strength 7 . 
The squares and circles are respectively the results of the 
FCI and CC calculations. The black dashed line is the an¬ 
alytic two-body result £2 for l-fl particles [57], and the red 
dotted line is the thermodynamic result £00 obtained from 
the GY-I-LDA approach [18]. The inset shows the difference 
£00 — £2 vs. 7 , which remains surprisingly small for every in¬ 
teraction strength. 


the presence of a strong transverse confinement. Along 
the axial direction the atoms are further confined by a 
harmonic potential of frequency w, and they interact by 
means of a short-range (contact) interaction of strength 
g. The corresponding Hamiltonian is then 
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In the absence of interactions, the total energy of the 
trapped gas is while the Fermi energy 

(defined as the energy of the first unoccupied level) is 
Ep = {N+l)fujj/2. The interaction strength may be suit¬ 
ably parametrized in terms of the dimensionless constant 
7 = {irg/y/N)/{hLoaz), where = sfhjrnw is the har¬ 
monic oscillator length. As the interactions are varied, 
the system continuously evolves from a Tonks-Girardeau 
(TG) gas of N strongly repulsive fermions (at 7 1) to 

a Lieb-Liniger (LL) gas of 7V/2 hard-core bosonic dimers 
(at 7 < - 1 ) [18]. 

Full configuration interaction (FCI) and CC 
methods. The exact solution of the many-body 
Schrodinger equation of a fermionic system can always 
be written as a linear combination of all Slater determi¬ 
nants that may be obtained from a complete set of one- 
particle basis functions. In a numerical solution, how¬ 
ever, only a finite number rif, of functions may generally 
be taken into account. The simplest approximation to 
the exact wave function is the Slater determinant |$) ob¬ 
tained by solving the mean-field Hartree-Fock (HF) equa¬ 
tions. The best possible solution in a given finite basis 


set can be obtained by expanding the wave function in all 
Slater determinants that may be obtained by replacing 
(i.e., exciting) one-particle functions in |$), and optimiz¬ 
ing variationally the coefficients in the resulting linear 
combination. This method is termed full configuration 
interaction (FCI) [58], or exact diagonalization, and it 
provides a strict upper bound to the exact ground-state 
energy. The computational cost of an FCI calculation 
scales as for nt ^ N [59], which forces a trade¬ 

off between the number of particles in the system and 
the number of basis functions needed to accurately 
describe it. At present, we are limited to « 50 for 
N = 6. One possible way to overcome this obstacle is to 
limit the number of excitations included in the FCI wave 
function. This leads to considerable savings of computer 
time, but unfortunately any truncated Cl calculation is 
size-inconsistent, in the sense that the total energy of two 
systems, not interacting with each other, is not guaran¬ 
teed to be the sum of the energies of the two systems. 

To overcome the size-consistency problem, the cou¬ 
pled cluster (CC) method was introduced, first in nuclear 
physics [60] and shortly after in quantum chemistry [45]. 
The CC wave function is given by the exponential Ansatz 

|vl/)=e^|$), (2) 

where T = Ti + T 2 + ... + Tn is the sum of all possible 
excitation operators T}^ replacing k HF one-particle or¬ 
bitals in the reference Slater determinant with k orbitals 
which are not present in |<I>). Due to the exponential form 
of the Ansatz (2), the CC method truncated to single 
and double excitations effectively includes various triply, 
quadruply, and higher excited determinants, since it con¬ 
tains, e.g., the products T 1 T 2 and T|. The CC method 
including all excitations is obviously equivalent to the 
FCI method, and its computational cost is equally high. 
However, truncated CC calculations are by construction 
size-consistent, and are much less time consuming than 
FCI with the same basis set size ni,. In the present work 
we use the coupled cluster method restricted to single, 
double, and non-iterative triple excitations, CCSD(T), 
whose computational complexity scales as [61]. The 
CCSD(T) method is very accurate for many properties of 
atoms and molecules, and it is now considered the golden 
standard of quantum chemistry [50]. 

We construct our many-body wave function by ex¬ 
panding it on the restricted basis containing the first rib 
single-particle eigenfunctions of the ID harmonic oscilla¬ 
tor, ())„(z) = H„(z/az) exp (-z^/ 2 a^)/y/ 2 ^bTla^^, with 
Hn{ ) an Hermite polynomial. A detailed study of the 
convergence with the type of excitations included in the 
coupled cluster model and the size of the basis set have 
been reported elsewhere [62]. For the purposes of this 
work, we simply note that the CCSD(T) method ap¬ 
pears to be the best choice for ID fermionic systems, 
in terms of the tradeoff between accuracy and computa¬ 
tional cost. As the energy is found to converge to the 
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Figure 2. Difference between the rescaled few-body energies 
and the analytical 2 -body result; filled circles (squares) are the 
CC (FCI) results, while the red dotted line is the GY-I-LDA 
result. The insets are vertical cuts across the main figure (i.e., 
at fixed 7 ), plotted as a function of 1/N. The bright green 
dots are the expected many-body limit, and the empty circles 
are the first order perturbationb approximation, Eq. (4). 


Figure 3. BCS pairing gap Av of a few-fermion system with 
N particles. Filled circles indicate CC results for increasing 
N, and the red squares are their extrapolation to A —>■ 00 . 
The red dotted line is the thermodynamic result, Eq. (5). 
Inset: chemical potential /rjv [in units of Nhuj] vs. 1/N, eval¬ 
uated at 7 = —Tx 12. 


exact value with a rate ^ l/.^/nb, all results presented 
here are obtained by extrapolating to the limit of infi¬ 
nite basis set with a quadratic interpolation vs. 
passing through the energies obtained for three values of 
rib up to 200. The FCI and CC calculations were per¬ 
formed with customized versions of the HECTOR [63] 
and ACESII codes [64], respectively. 

Results. We start by considering the ground state en¬ 
ergy Ejq of a balanced ensemble of N interacting fermions 
in the trap. In order to appropriately compare ensembles 
with different numbers of particles, in Fig. 1 we plot the 
dimensionless quantity as a function of 

the rescaled interaction strength 7 . The result is analytic 
for the simplest case of two (1 -I- 1) particles, and is the 
solution of the implicit equation [57] 

1 ^ F(3/2 - E 2 / 2 ) 

g 72(^2 -i)f(i-f; 2 / 2 )' 

In the thermodynamic limit of an infinite number of par¬ 
ticles, instead, the analytic result Ego can be obtained 
by applying the local density approximation (EDA) to 
the solution of the Gaudin-Yang (GY) integral equations 
describing a homogeneous gas [18]. We start by noticing 
that the results in the two extreme limits are actually 
surprisingly close to each other. As shown in the inset 
of Fig. 1, the many-body and two-body rescaled ener¬ 
gies differ by less than 0.05 over the complete range of 
interaction strengths. This poses a serious challenge to 
the numerical calculation of the energies. Nonetheless, 
we see that the results of both CC and FCI lie nicely be¬ 
tween the two curves, thereby on one side showing their 
accuracy, and on the other providing a further confirma¬ 
tion of the validity of EDA for computing the energy of 


a trapped ID gas. 

The size of the computational space diverges exponen¬ 
tially with the number of particles within exact diagonal- 
ization, so that with this method we could obtain con¬ 
verged results only for very small systems, containing at 
most three pairs of atoms {N = 6 ), compatibly with what 
was recently published in Ref. [35]. On the other hand, 
only a carefully chosen subset of the total Hilbert space 
is retained in the CC calculations, so that this method 
allows us to investigate much larger systems, even up to 
40-1-40 particles. The range of interaction strengths we 
are able to explore however slowly shrinks with the size of 
the system N, as the complexity of the calculation scales 
with the bare interaction strength g, rather than with its 
rescaled counterpart 7 oc g/ y/N. 

To investigate in detail the continuous crossover from 
few- to many-body systems, in Fig. 2 we show our results 
after subtracting the analytical two-body contribution 
82 - Even after the subtraction, the energies are shown 
to converge remarkably fast to the thermodynamic limit 
(red dotted line). While our CC and FCI results closely 
match for 2 - 1-2 particles, at this level of precision one no¬ 
tices that the FCI results for 3-1-3 particles start to devi¬ 
ate from the expected behavior: as an example, for 7^2 
the 3-1-3 FCI results (orange squares) unphysically cross 
the 2 - 1-2 results (blue symbols), even after a very time- 
consuming calculation (a week of CPU time to get one 
not-yet-converged FCI point, compared to two hours for 
a converged CCSD(T) point). The insets of Fig. 2 show 
the rescaled CC few-body energies plotted at fixed inter¬ 
action strengths versus 1/iV. For large N the few-body 
results smoothly extrapolate to the GY-I-EDA thermo¬ 
dynamic limit (green dots), while for small g (i.e., small 
7 and small N) the results match the first order pertur- 
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bation result 

/ OO 

dznl + 0{g^), (4) 

-OO 

where no = is the density of a gas of 

N/2 identical fermions. We note in passing that, in the 
limit of small g and large N, the latter equation recovers 
the expected weak-coupling LDA result, En = + 

4gN^/y{57T^a,) + Oig^) [ 22 ], 

We turn now to consider the BCS pairing gap, which 
for this few-body system we define as Ajv = E]\f — 
{En+i + En-i)/^, for odd values of N (with = 
Ni + 1). The results of our CC calculations are shown 
in Fig. 3. The BCS pairing gap equals half the spin gap, 
and in the thermodynamic limit for a homogeneous sys¬ 
tem it is identically zero for repulsive interactions, while 
it has a characteristic non-analytic behavior for small at¬ 
tractive interactions. The GY result for a homogeneous 
gas [ 20 ] may be adapted to describe a trapped system 
by replacing the homogeneous Fermi momentum Tm/2 
(n being the total density) with its value at the center 
of the trap, kp ^ y/N /qz (for N ^ 1). The specific 
choice of the momentum at the center of the trap can 
be justified by the fact that, even if the extra particles 
are added at the edges of the cloud, their presence will 
cause an overall reorganization of the gas density profile, 
so that the resulting energy gap will be sensitive to the 
typical momentum kp. This procedure yields 

(y) (5) 

Our few-body results are, as expected, analytic across 
the non-interacting point. However, one can clearly see 
how the progressive build-up of the Fermi sea gives rise 
to the expected non-analytic behavior in the weak in¬ 
teraction limit, as predicted by Eq. (5). The inset of 
Fig. 3 shows instead the behavior of the chemical po¬ 
tential /rjv = Ajv+i — Ejv, which is characterized by a 
pronounced even/odd effect due to the shell structure of 
the trap. At odds with the energy, the approach to the 
thermodynamic behavior is actually very slow for both 
the pairing gap and chemical potential. 

Finally, of great interest is the density profile of the 
trapped cloud. The density at a point zq is obtained 
as the expectation value of the operator ~ ^o) 

within the FCI calculations, and in the CC approach 
from the Hellmann-Feynman theorem by adding a per¬ 
turbation of the form X(j3i(zo)(f>j{zo) to the ij-th element 
of the one-particle Hamiltonian matrix, and taking the 
derivative of the energy with respect to A. The density 
profiles we obtain are shown in Fig. 4. These become 
broader as the interaction strength grows increasingly 
more repulsive, and display a regular series of peaks. 
One finds exactly N/2 peaks in the LL limit of strong 
attraction (7 ^ — 1 ), where pairs of fermions be¬ 
come tightly-bound bosonic hard-core dimers, and the 
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Figure 4. Density of a balanced two-component Fermi gas. 
The solid lines in the left panel display FCI results for 3-1-3 
particles, while the ones in the right panel show CC results 
for 15-1-15 atoms. The dashed lines are the GY-I-LDA results, 
and the black dotted lines are analytic results for 7^—1 
(LL), 7 = 0 (0), and 7 » 1 (TG). 

density approaches till = ^ with ^i(z) the 

wavefunction of the i-th level for a particle of mass 2 m. 
N/2 peaks are present also in the non-interacting case, 
where two distinguishable fermions occupy the same or¬ 
bital, and the density becomes 2no. The number of peaks 
then smoothly evolves to N in the TG limit (7 ^ 1) 
of “fermionized” fermions, where due to the strong re¬ 
pulsion even distinguishable fermions must occupy dif¬ 
ferent one-particle levels, and the density approaches 
^TG = In fhs thermodynamic limit, the 

GY-I-LDA analysis [18] predicts a density with the typi¬ 
cal Thomas-Fermi (TF) profile of an inverted parabola, 
np-piz) cx (1 — z^/i?^p), where the radius Rtf varies 
between s/N/^az in the LL limit and \/2Naz in the 
TG limit, being exactly equal to '/NcLz for an ideal 
gas. The two approaches used here prove to be some¬ 
how complementary. The FGI method allows us to con¬ 
sider very strong repulsive interactions and enter the ex¬ 
treme “fermionized” regimes: the hard-core bosonic LL 
gas, and the fermionized atomic TG gas. Such strong 
repulsion is out of the reach of present CG calculations, 
since in this regime both the Hartree-Fock and the corre¬ 
lation energies diverge. However, with CC we are able to 
address much larger systems and explore beyond mean- 
field corrections. In particular, this method gives direct 
access to non-perturbative effects such as the progressive 
build up of a non-analytic behavior for the pairing gap. 

Summarizing, we presented a detailed analysis of the 
static properties of a two component Fermi gas in a ID 
harmonic trap. We computed with high accuracy the en¬ 
ergy, the chemical potential, the pairing gap, and density 
profiles for N ranging from a few to a few tens, and for 
a broad range of interaction strengths, well beyond the 
mean-field regime. Our predictions for £jv, A*tv and A a? 
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may be tested directly by RF spectroscopy or by tun¬ 
neling measurements, as done in Refs. [13, 14]. The CC 
method proved to be extremely well suited to studying 
harmonically trapped fermions, and we foresee that it 
will equally well describe more complex potentials, such 
as double wells or microtrap arrays [15], and dipolar sys¬ 
tems. The method could moreover be generalized to 2D 
and 3D, or to few strongly-interacting bosonic atoms. 
These studies would, in particular, be crucial to under¬ 
stand whether the rapid convergence of the system’s en¬ 
ergy found here is mainly due to the short-range charac¬ 
ter of the interactions, to the ID confinement, or to the 
harmonic external potential. 
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